Butterflies Solve Bidiagonal Toeplitz Systems

Here and hereafter T = (ti,j) n−1 i,j=0, ti,j = 1 for i − j = 0, ti,j = c 6= 0 for i − j = 1, and ti,j = 0 otherwise, x = (xi) n−1 i=0 , and b = (bi) n−1 i=0 . Note that the system is scaled so that the main diagonal is composed exclusively of ones, with no loss of generality. These systems are at the heart of problems as diverse as cubic spline and Bspline curve fitting [3], [11], preconditioning for iterative linear solvers [2], [13], computation of photon statistics in lasers [8], computational fluid dynamics [26], solution of multidimensional diffusion computations [9], solution of neuron models by domain decomposition [14], and more. We propose a parallel algorithm based on the fast Fourier transform (FFT) to solve such systems. Cooley and Tukey [4] introduced the FFT in 1965 as a mechanism for fast computation of the discrete Fourier transform (DFT) on computers. General parallel architectures tend to suffer from high latency and restrictive bandwidth for communication between processing units and/or clusters of processing units and dedicate effort to fetching and decoding instructions. Specialized hardware to perform the FFT in parallel alleviates these impediments to efficient parallel computation. A simple set of computations at the heart of the FFT can be carried out by a circuit known as the Butterfly. A Butterfly circuit as depicted in Figure 1 accepts two inputs. One input is multiplied by what is referred to as a twiddle factor, which are nth roots of unity for an FFT of dimension n. Then in parallel this product is both added to and subtracted from the other input. The two resulting values are the output of the Butterfly circuit. FFT processors generally incorporate log2 n stages of n/2 parallel Butterfly circuits to provide pipelined computation of an n point DFT.